/*  sdfrcbrt.cpp  by K.Tsuru  */
// function ID 3011     DRADIX only
/*****************************************************************
SDouble class
reciprocal cubic root of x
[Algorithm]
Get a solusion of an equation 1/(x^3) - a = 0 by the
Newton's method with "doubling effective figures"(DEF).
See Sqrt() for detail.
Let d(n) = {x(n) - a*x(n)^4}/3 ... faster than x(n){1-a*x(n)^3}/3
    x(n+1) = x(n) + d(n)
******************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
static const char* const func = "RecCbrt";
SDouble RecCbrt(const SDouble& a){
    if(a.Type() != a.REAL) a.SetError(a.RADIX_ERR, func, 3011);
    if(a.Sign() == 0) a.SetError(a.DIVIDED_BY_ZERO, func, 3011);
    RealSize C;
    uint max_sz = a.MaxSize();
    SDouble w(Dabs(a)), one(1.0);
    double w0;

    if(a.IsOne()){  // |a| = 1.0
        if(a.Sign() < 0) one.ChangeSign();
        return one;
    }
    int e = w.NetRdxExp(), q, r, pow10 = 0;
// a = d*R^e = (d*R^q)*R^(3*r) , e = q + 3r
    r = e/3;
    q = 3*r;
    if(e - q > 0){ r++; q += 3; }
    w.MultPowRdx(-q);   // Here 1/R^3 <= w < 1.0.
/*set an initial value y0*/
    w0 = doubleD(w);
    while(w0 < 1.0e-3){ // 0.001 <= w0 < 1.0
        w0 *= 1000.0;   pow10++;
    }
    if(pow10) w.MultPow10(3*pow10);
//enter the irreducible size mode
    SDouble y(a.REAL, max_sz), delta(a.REAL, max_sz), z(a.REAL, max_sz);
    y = pow(w0, -1.0/3.0);  // y has about DOUBLE_FIG figures.
/***************************************************************/
    int c = 0, itrmax = howpow2( (DFIGURES*max_sz)/DOUBLE_FIG+1 ) +6;
    uint ef = (DOUBLE_FIG*2u)/DFIGURES, fig = C.EffFigures();
    int fullPrec = 0;
    delta = w0 - w;
    if( delta.Sign() ){
    // If w0 is very close to double value, higher precision is necessary.
        ef = max( ef, 2u*(uint)abs( delta.NetRdxExp()) );
    }
//enter fiexd point mode
    y.FixedPoint(y.RdxExp());
    if(ef > fig) ef = fig;
    do{
        if( (ef = C.SetEffFig(ef)) >= fig) fullPrec = 1;
        z = y*y;    z = z*z;         // z = y^4
        delta = y - w*z;
        delta = DsDiv(delta, 3);    // delta /= 3
        y += delta;
        c++;
        ef *= 2;
        if(ef >= fig) ef = fig;
        C.SetEffFig(0);
    } while( ( !delta.IsHidden() && (c < itrmax) ) || !fullPrec );
    y.iterationCount = c;
#ifndef NDEBUG
    if(c == itrmax) y.SetError(y.FATAL, func, 3011);
#endif
    y.PointFree();
    y.Reform(3011);
/*******************************************************************
[Reference]
A revision may not be necessary.
Let d = 1/(y'^3) - w.
If d != 0
    y = w^(-1/3) = {1/(y'^3) -d}^(-1/3) ~= y'{1 + d*(y'^3)}/3
    = y' + (y'^4)*d/3.
********************************************************************/
#if 0
    if(y.PreferSpeed() == OFF){
        delta = Dpow(y, -3) - w;    //Almost 'delta' has only one figure.
        if(delta.Sign()){
            z = y*y;    z = z*z;    // z=y^4
            delta = delta*z;
            y += DsDiv(delta, 3);
        }
    }
#endif

    if(y.Verify()){
        z = y*y;    z *= y; // z =y^3
        delta = one - z*w;  // delta = 1 - w*y^3
        if(!delta.IsMLT(w)){
            y.SetError(y.VERIFY, func, 3011);
        }
    }
    if(r) y.MultPowRdx(-r);
    if(pow10) y.MultPow10(pow10);
    if(a.Sign() < 0) y.ChangeSign();
    return y;
}